##########################################################################################

library(data.table)
library(optparse)
library(ArchR)
library(BSgenome.Hsapiens.UCSC.hg38)

##########################################################################################
option_list <- list(
    make_option(c("--comine_data_file"), type = "character"),
    make_option(c("--peak_gene_file"), type = "character"),
    make_option(c("--cluster_type_file"), type = "character"),
    make_option(c("--cluster"), type = "character"),
    make_option(c("--out_path"), type = "character") 
)

if(1!=1){
    
    ## 整合atac和rna的文件
    comine_data_file <- "~/20231121_singleMuti/results/qc_atac_v3/somatic/testis_combined_peak.combineRNA.qc.Rdata"

    ## 细胞类型注释文件
    cluster_type_file <- "~/20231121_singleMuti/config/cluster_celltype.csv"

    cluster <- "cluster0"

    ## peak-gene文件
    peak_gene_file <- "~/20231121_singleMuti/results/celltype_plot/peak2gene/somatic/peakToGeneHeatmap_LabelClust_k25.tsv"

    ## 输出
    out_path <- "~/20231121_singleMuti/results/subcell"

}

###########################################################################################
parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

comine_data_file <- opt$comine_data_file
peak_gene_file <- opt$peak_gene_file
cluster_type_file <- opt$cluster_type_file
cluster <- opt$cluster
out_path <- opt$out_path

dir.create(out_path , recursive = T)

###########################################################################################

a <- load(comine_data_file)
atac_proj <- testis_combined_peak_combineRNA

dat_peak_gene <- data.frame(fread(peak_gene_file))

## 分细胞类型跑
cluster_type <- fread(cluster_type_file , header = F)
celltype_name <- subset( cluster_type , V1 == cluster )$V2

###########################################################################################
## 提取细胞类型
subCells <- getCellNames(testis_combined_peak_combineRNA)[as.character(testis_combined_peak_combineRNA@cellColData[["cell_type"]]) %in% celltype_name]

## 提取特定细胞的
sub_proj <- subsetArchRProject(
  ArchRProj = testis_combined_peak_combineRNA,
  cells = subCells,
  outputDirectory = out_path,
  dropCells = TRUE,
  force = TRUE
)

###########################################################################################
## 基于chromvar鉴定motif
projHeme5 <- sub_proj

#projHeme5 <- addMotifAnnotations(ArchRProj = projHeme5, motifSet = "cisbp", name = "Motif" , force = TRUE)

#projHeme5 <- addBgdPeaks(projHeme5 , force = TRUE)

#projHeme5 <- addDeviationsMatrix(
#  ArchRProj = projHeme5, 
#  peakAnnotation = "Motif",
#  force = TRUE
#)

###########################################################################################
## scATAC-seq 和 scRNA-seq 数据识别峰到基因链接
#projHeme5 <- addPeak2GeneLinks(ArchRProj = projHeme5,useMatrix = "GeneExpressionMatrix")

## 提取感兴趣的peak集合
tmp_peak <- data.frame(projHeme5@peakSet)
tmp_peak <- paste0(tmp_peak$seqnames , ":" , tmp_peak$start , "-" , tmp_peak$end)
use_peak <- tmp_peak[tmp_peak %in% unique(dat_peak_gene$peak)]
tmp_peak_index <- which(tmp_peak %in% use_peak)

## peak只用存在peak-gene关系的peak
projHeme5@peakSet <- projHeme5@peakSet[tmp_peak_index,]
## 使用自己定义的peak集合
projHeme5 <- addPeakSet(projHeme5,projHeme5@peakSet,force = TRUE)
projHeme5 <- addPeakMatrix(projHeme5)

projHeme5 <- addMotifAnnotations(ArchRProj = projHeme5, motifSet = "Cisbp" , name = "Motif",force = TRUE)
projHeme5 <- addBgdPeaks(projHeme5 , force = TRUE)
projHeme5 <- addDeviationsMatrix(
  ArchRProj = projHeme5, 
  peakAnnotation = "Motif",
  force = TRUE
)

atac_proj <- projHeme5


###########################################################################################
out_file <- paste0( out_path , "/" , cluster , ".combineRNA.motif_peak2gene.Rdata" ) 
save( atac_proj, file = out_file )
saveArchRProject( atac_proj , out_path )
